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Abstract. We investigate the connection between the linear harmonic oscillator 
equation and some classes of second order nonlinear ordinary differential equations of 
Lienard and generalized Lienard type, which physically describe important oscillator 
systems. By using a method inspired by quantum mechanics, and which consist 
on the deformation of the phase space coordinates of the harmonic oscillator, we 
generalize the equation of motion of the classical linear harmonic oscillator to several 
classes of strongly non-linear differential equations. The first integrals, and a number 
of exact solutions of the corresponding equations are explicitly obtained. The devised 
method can be further generalized to derive explicit general solutions of nonlinear 
second order differential equations unrelated to the harmonic oscillator. Applica¬ 
tions of the obtained results for the study of the travelling wave solutions of the 
reaction-convection-diffusion equations, and of the large amplitude free vibrations 
of a uniform cantilever beam are also presented. 
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1. Introduction 


An important equation in mathematical physics is the Lienard type 
ordinary second order nonlinear differential equation of the form [1] 

x(t) +f(x)x(t)+g(x) = 0, (1) 

where a dot represents the derivative with respect to the time t, and / 
and g are arbitrary C 1 functions of x. Together with its generalization, 

* email: t.harko@ucl.ac.uk 
t email: stslsd@mail.sysu.edu.cn 



© 2015 Kluwer Academic Publishers. Printed in the Netherlands. 


Lienard_JENG.tex; 30/06/2015; 1:39; p.l 



2 


T. Harko and S. D. Liang 


the Levinson-Smith type equation [2] 

x{t) + / (x,x) x(t) + g{x) = 0, (2) 

where / is an arbitrary function of x and x, the Lienard equation 
plays an important role in many areas of physics, biology and engi¬ 
neering [3]. The mathematical properties of these types of equations 
have been intensively investigated from both mathematical and physi¬ 
cal point of view, and their study remains an active field of research in 
mathematical physics [4] - [10]. 

One of the major scientific and engineering applications of the gen¬ 
eralized Lienard type equations is the mathematical study of a large 
number of non-linear oscillations, which can be described by 

x + f (x) g(x) + h(x) = 0, (3) 

where /, g and h are arbitrary C 1 functions. In many practical situa¬ 
tions the function / (x) can be approximated as 

f (x) = ax 3 + fix 2 + jx + 5, (4) 

where a, /3, 7 , 5 are constants [11, 12]. From a physical point of view the 
Lienard equation can be interpreted as the mathematical generalization 
of the equation of the damped oscillations, 


x + 'jx + uj 2 x = 0, (5) 

with 7 = constant and u 2 = constant, respectively [13]. For 7 = 0 
we obtain the equation of the linear harmonic oscillator, which repre¬ 
sents one of the fundamental equations of both classical and quantum 
physics. Generally, a linear oscillation can be described by the second 
order ordinary differential equation [13] 

x + f(t)x + g(t)x = 0 . ( 6 ) 

The second order linear equation describing the oscillations of a com¬ 
plex physical system can be reduced to a Riccati equation [14, 15, 16]. 
Hence the solutions of the Riccati equation corresponding to the given 
oscillator can be used to study the linear oscillations of classical sys¬ 
tems. By using N. Euler’s theorem on the integrability of the general 
anharmonic oscillator equation [17], three distinct classes of general 
solutions of the highly nonlinear second order ordinary differential 
equation 

(ft nr r] nr 

^ + Mt)- + f 2 (t)x + f 3 (t)x n = 0, (7) 
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were presented in [18]. A class of exact solutions for the Lienard type 
ordinary non-linear differential equation was obtained in [19]. In or¬ 
der to obtain the solutions, the second order Lienard type differential 
equation is transformed into a second kind Abel type first order dif¬ 
ferential equation. With the use of an exact integrability condition 
for the Abel equation (the Chiellini lemma) [20], the exact general 
solution of the Abel equation can be obtained, thus leading to a class 
of exact solutions of the Lienard equation, expressed in a parametric 
form. The Chiellini integrability condition was also extended to the case 
of the generalized Abel equation. As an application of the integrability 
condition the exact solutions of some particular Lienard type equa¬ 
tions, including a generalized van der Pol type equation, were explicitly 
obtained. The connection between dissipative nonlinear second order 
differential equations and the Abel equations, which in its first kind 
form have only cubic and quadratic terms, were considered in [21]. By 
using the integrability criterion of Chiellini, the general solution of the 
corresponding integrable dissipative equations was obtained. Through 
a modified factorization method nonsingular parametric oscillators, 
which are Darboux related to the classical harmonic oscillator, and 
have periodic dissipative/gain features, have been identified [22], Note 
that the same method was also applied to the upside-down (hyperbolic) 
” oscillator”. For these oscillators the obtained Darboux partners show 
transient underdamped features. Some basic results in the integrability 
of the Abel equations were reformulated and expanded in [23]. 

An interesting nonlinear oscillator, which is described by a Lienard 
type equation, and which has a number of interesting characteristics, 
is 


k 2 o 

x + kxx + —x + Aix = 0, 
9 


k. Ai = constant. 


( 8 ) 


The mathematical properties of this equation have been studied ex¬ 
tensively in [24]-[30]. Eq. (8) can also be considered as a generalized 
Emden-type equation. It has the interesting property that for Ai > 0 
nonisolated periodic orbits can be explicitly obtained. These periodic 
orbits, which are of conservative Hamiltonian type, have the interesting 
property that the frequency of the damped oscillations is completely 
independent of their amplitude. Moreover, they continue to have the 
same value as the one corresponding to the linear harmonic oscillator 
[24], In [25] it was shown that Eq. (8) is integrable for Ai = 0. Eq. (8) 
can be solved either explicitly, or by quadratures. It also turns out that 
it can be equivalently described in terms of some time-independent 
Hamiltonian functions. A particular case of Eq. (8), 


x + 3xx + x 3 + lo 2 x = 0, 


(9) 
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was studied in [26]. Eq. (9) has also the uncommon property that the 
frequency of oscillation of the corresponding oscillator is independent 
of its amplitude. This property is specific for linear harmonic oscil¬ 
lators, and it is not a general characteristic of nonlinear oscillator 
equations. An oscillator possessing the property of independence of 
frequency from amplitude is also known as an isochronous oscillator 
[27]. The isochronicity property of Eq. (9) can be related to its main 
characteristic, namely that it can be transformed through a nonlocal 
transformation [26] 

U(t) = x(t)efi x W dt ',. (10) 

to the linear harmonic oscillator equation, 

U(t) +uj 2 U(t) = 0. (11) 

A more general non-local transformation, 

U(t) = x(t)eti f W» dt ’ , ( 12 ) 

as applied to the linear harmonic oscillation equation, generates a non¬ 
linear Lienard type differential equation of the form [26] 

x + (2/ + xf) x + (/ 2 + to 2 ) x = 0, (13) 

where the prime denotes differentiation with respect to x. One can also 
consider more general nonlocal transformations, defined as [26] 

U(t) = g (x(t)) efi , (14) 

leading to the equation 

g" 2 9f 9 2 f 2 9 

X + x 2 + ^—x + 2 fx + —u? + J —r- = 0, (15) 

9 9 9 9' 

where g' = dg/dx. A simple and efficient procedure to construct higher 
dimensional isochronous Hamiltonian systems, called the O-modified 
procedure, was developed in [28]. By using this procedure it was shown 
that a class of singular Hamiltonian systems is equivalent to constrained 
Newtonian systems. The procedure was first proven for two dimensional 
dynamical systems, and then extended to the case of N-dimensional 
systems. It was also shown that, depending on the choice of the sys¬ 
tem parameters, the considered systems admit two types of solutions: 
isochronous, and amplitude independent quasiperiodic solutions, re¬ 
spectively. The quantization of the nonlinear oscillators of type (5) was 
considered in [29] and [30], respectively. 
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A very powerful method of solving the differential equation, 

y"{x) + f 3 (x, y)y 2 + f 2 (x,y)y' + fi(x,y ) = 0 , (16) 

where fi(x,y), i = 1,2,3 are arbitrary functions, was proposed in [31]. 
This method is as follows: assume that / 3 = f 3 (y) and fi = f- 2 (x). 
Then, provided that has the form 


h(x,y) 


q(x)Z(y) - exp 



f 2 (x)dx 


C(y ), 


(17) 


where q(x) is arbitrary, and Z’(y) := l~f 3 Z, C\y) := - (/ 3 + 3/Z) C(y), 
the general solution of Eq. (16) is y = J(r, s), where r and s are variable 
independent solutions of the linear equation 

y”{x) + f 2 (x)y' + q{x)y = 0 , (18) 

and the function J is any solution of the system 

J rr = S 2 C(J) - h(J)J 2 ., J rs = ~rsC(J) - f 3 (J)J r J s , (19) 

J ss = r 2 C(J)- h(J)Jl J r = [Z(J)-sJ s }/r. (20) 

In [32] it was shown that Eq. 16) can be transformed into a first order 
equation, instead to a second order linear equations. Using this trans¬ 
formation, a first integral of Eq. (16) can be found. The results of [32] 
were generalized in [33], where a method for finding the complete first 
integral of the Eq.(16) was obtained, in the case when the coefficients 
fi, i = 1, 2, 3 fulfil a certain condition. 

A method for obtaining a first integral of a Lienard type equation 
of the form 

x + — \xf(x)] x + u 2 x + ^ — = 0, (21) 

x ax x 

where f(x) is an arbitrary C 1 function, was introduced in [34]. The 
method is based on the introduction of a first order linear ’’generating 
equation” for a complex function X = X(x,x), which is equivalent 
to the equation of the linear harmonic oscillator, and by considering 
a ’’deformation” of this equation by means of the substitution x —> 
x + f(x). 

It is the purpose of the present paper to generalize the approach 
introduced in [34], First, by introducing, in analogy to the quantum 
mechanical treatment, two ’’creation” and ’’annihilation” functions a = 
a(x,x,u) and a + = a + (x,x,lu), lo = constant, respectively, we define 
a generating function X = X (a, a + ), satisfying a given particular dif¬ 
ferential equation. For a specific choice of a and a + it follows that x 
satisfies the simple linear harmonic oscillator equation, whose solution 
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can be therefore obtained by using only the properties of the functions 
a and a + . By introducing a general ’’deformation” of the space-phase 
coordinates (x, x ) so that 

x —> x = x + g (f , x, x), (22) 

and 

x x = x + f (t, x, x) , (23) 

with / and g being arbitrary C 1 functions of f, x, and x, we construct 
the corresponding deformations of the functions a (x, x, u) — >• a (x, x, w) 
and a + (x,x,u ;) —> a + (x, x, w), respectively. 

By imposing the condition that X = X (a, a + ) satisfies the same 
equation as in the case of the harmonic oscillator, a large class of 
non-linear Lienard, or generalized Lienard, type differential equations 
can be generated. These equations represent a non-linear extension of 
the harmonic oscillator equation, and their first integral can be easily 
obtained. Therefore the exact solution of several distinct classes of 
ordinary second order non-linear differential equations can be obtained 
in an exact analytical form. 

The present paper is organized as follows. In Section 2 we review the 
problem of the linear harmonic oscillator, by using methods from quan¬ 
tum mechanics, and applying them to the classical case. The solution 
generating first order differential equation, and its solution are explic¬ 
itly presented. In Section 3 we consider the phase space deformation of 
the coordinates x and x. The resulting differential equations, and their 
first integrals are explicitly presented. The exact solution of some non¬ 
linear differential equations are also obtained. We present a number 
of possibilities of further generalizing our results in Section 4. Some 
physical and engineering applications are briefly analyzed in Section 5. 
We discuss and conclude our results in Section 6 . 


2. The linear harmonic oscillator 

The equation of motion of the linear harmonic oscillator is given by 

x(t) + ui 2 x(t ) = 0, cu = constant. (24) 

Based on the analogy with quantum mechanics [35] we introduce two 
complex ’’creation” and ’’annihilation” functions a(t ) and a + (t), defined 
as 

a(t ) := x(t) — i ujx(t), (25) 

and 

a + (t ) = a*(t ) := x(t) + it ox(t) (26) 
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respectively, where * denotes complex conjugation, and x{t ) is a solu¬ 
tion of Eq. (24). Due to their definition, and the equation of motion of 
the harmonic oscillator, the functions a and a + have the properties 

H(t ) = a(t)a + (t) = x 2 (t) + u 2 x(t) 2 , (27) 


and 

dH ^ = ^ [a(t)a+(t)] = a(t)a + (t) + a{t)a + {t) = 0, (28) 

respectively. We introduce now the solution generating function X 
defined as 

a(t ) x(t ) — i oox(t) 


X(t) = 


(29) 


a + (t ) x(t) + iojx(t) 

It follows that Eq. (24) is equivalent with the solution generating 
equation [34], 

X(t) = -2i u>X(t). 

Indeed, Eq. (30) gives 

aa + — aa + = —2icuaa + , 
or, with the use of Eq. (28), 


a = —lcua, 


(30) 

(31) 

(32) 


from which Eq. (24) can be immediately obtained. 

Eq. (30) has the general solution 

X(t) = e ~ 2 ^ t+a \ 

where a is an arbitrary constant of integration. Therefore from 
a(t) x(t) — iujx{t) 


X(t) = 


_ g —2i (ut+a) 


a + (t) x(t) + iux(t) 
we obtain the first order linear differential equation 


(33) 


(34) 


x(t) = lo cot (uot + a) x(t), (35) 

with the general solution x(t) = A sin (c ut + a), where A is an arbitrary 
constant of integration. Therefore we have obtained the general solution 
of the equation of motion of the linear harmonic oscillator by using the 
solution generating Eq. (30), and the definitions of the ’’creation” and 
’’annihilation” functions a(t) and a + (f), respectively. In [36] and [37] an 
alternative factorization of the Hamiltonian of the quantum harmonic 
oscillator was considered. This factorization leads to the construction of 
a specific two-paranreter self-adjoint operator, from which the standard 
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harmonic oscillator equation can be recovered in certain limits of the 
model parameters. 

In the case x = 0, x ^ 0, we have X(t) = 1, a = a + . Then, from 
Eq. (33) it follows that the oscillating system reaches the position x = 0 
at time moments given by 

t n = n <E N. (36) 

OJ 


3. Generalizing the linear harmonic oscillator equation by 
deforming the function a(t ) 

In the present Section we will generalize the form of the functions 
a and a + , and consequently of the solution generating function X , 
by ’’deforming” the phase - space coordinates x and x according to 
Eqs. (22) and (23). The deformation of the phase space coordinates 
will induce a deformation of the ” creation” and ” annihilation” functions 
a (x, x) and a + (x, x) of the linear harmonic oscillator. Therefore under 
the deformation of the coordinates the ” creation” function a transforms 
as 

a (. x , x,uj) a (x + g (f, x, x ), x + / (t, x, x ), u). (37) 

We will continue to assume that the ’’deformed” operator X, given 
by 

^ = _a_ = x(t) + / (t, x, x) - itu \x(t) + g (t, x, .7)] 
a + x{t) + f (t, x, x) — iuj [x(t) + g (t, x, x)\ ’ 

satisfies the solution generating equation Eq. (30) of the simple har¬ 
monic oscillator, 

X + 2i ujX = 0, (39) 

with the general solution given by 

X(t) = e ~ 2i(u}t+a) . (40) 

The equation 

d = d + e~ 2i(u}t+a \ (41) 

gives the first integral of Eq. (39). Let F (t,x,x,x) = 0 a second order 
ordinary differential equation. A relation of the form 4> ( t, x, x) = C, 
where C is a constant, and x is a solution of the equation F (t, x, x, x) = 
0, is called a first integral of F ( t,x,x,x ) = 0 [38]. 

The procedure based on the ’’creation”, ’’annihilation” and gener¬ 
ating functions allows us to obtain the first integrals, and therefore 
the exact solutions of a large class of ordinary non-linear differential 
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equations, including Lienard and generalized Lienard type equations. 
The most general deformation of the function a(t) involves the pres¬ 
ence of two arbitrary C 1 functions / and g. The case / (t, x, x) = 0, 
g (t, x, x) = 0 corresponds to the case of the simple harmonic oscilla¬ 
tor. For g(t,x,x ) = 0, and f(t,x,x) = x we reobtain the equation 
presented in [34], 

If the function g ( t,x,x ) satisfies the conditions 

x + g (t, x, x) = 0, x + / ( t , x, x) / 0, (42) 

similarly to the case of the simple harmonic oscillator it follows that 
X = 1, and the non-linear system passes through the point x+g (t, x, x) = 
0 at the same time intervals as the simple harmonic oscillator passes 
through the point x = 0 . 

In the following we will assume that / and g are independent C l 
functions of the independent coordinate f, of the function x, and of 
x, respectively. By adopting specific forms for the functions / and g 
we will investigate the corresponding classes of equations, and their 
solutions. 


3.1. a(t) = x(t) + f(t) — i oj [ x(t) + g(t)j 


As a first case for the generalization of the equation of the linear har¬ 
monic oscillations we consider the equation obtained by deforming the 
phase space coordinates x(t) and x(t ) so that x(t) —>• x{t) + g(t) and 
x(t) —>• x(t) + f(t), where f(t) and g(t) are arbitrary C 1 functions of 
time. Therefore the function a(t) becomes 


a(t) = x(t) + f{t) - ico [x(t) + g(t)} . (43) 


Then the solution generating equation Eq. (30) is equivalent with the 
second order non-linear differential equation 


[x(t) + g(t)\ x(t) + [f(t) - g(t)\ x(t) + J 2 [x(t) + g(t)f + 
f(t) [x(t) + g{t)] + [f{t ) - g(t)} f{t) = 0. (44) 

Eq. (44) has the first integral, 


x = oj cot (uit + a ) [x(t) + g(t)] — f(t ), 

which follows from X(t) = a(t)/a + (t ) = e ~ 2l (^*+«). 
Therefore the general solution of Eq. (44) is given by 


x(t) 


sin (cut + a) 


A + 


a; cot (ut + a)g(t) - f(t ) ' 

sin (cot + a) 


(45) 


(46) 
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where A is an arbitrary constant of integration. 

3.1.1. Particular case 1: g(t ) = 0 
If g(t ) = 0, Eq. (44) takes the form 

x(t)x(t) + f(t)x(t ) + ai 2 x 2 (t) + f(t)x(t) + f 2 (t) = 0, (47) 

with the first integral given by 

x(t) = ai cot (ait + a) — f(t). (48) 

Therefore the general solution of Eq. (47) is given by 

x(t) = sin (ait + a) A — f - ^ - -dt . (49) 

w v _ J sin (ait + a) y ’ 

As a simple application of the solution given by Eq. (49) we consider 
the case f(t) = fo sin (ait + a) , with /q constant. Then Eq. (47) takes 
then the form 

x(t)x(t) + fo sin (ait + a) x(t) + ai 2 x 2 (t) + foai cos (ait + a) x(t) + 


fo sin 2 (ait + a) = 0, (50) 

with the general solution given by 

x(t) = (A — fot) sin (ait + a). (51) 

3.1.2. Particular case 2: f(t ) = 0 
If /(f) = 0, Eq. (44) reduces to 

[x(t) + s(f)] x - g(t)x(t ) + a? [x(t) + g(t )] 2 = 0, (52) 

with the general solution given by 

, , . . . I" . f 01 cos (ait + a) . . , 1 . 

x(t) = sm(o;t + a) A+ -- g(t)dt . (53) 

|_ J sin (ait + a) 

If the function g(t) has the simple functional form g(t ) = go sin n (ait + a), 
go = constant, n € N > 1, Eq. (52) takes the form 

[x(t) + go sin” (ait + a)] x(t) — aingo sin ” -1 (ait + a ) cos (ait + a) x(t) + 
ai 2 [; x(t) + go sin" (ait + a )] 2 = 0, n € N > 1. (54) 

The general solution of Eq. (54) is given by 


/x . , . sm 1 (ait + a) _ 

x(t) = sm(ait + a) A + go --- , n € N > 1. (55) 

n — 1 
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3.2. a(t) = x + f(x) — vjj [.t + g{x)\ 

As a second case of the deformation of the function a we consider that 
a(t) is of the form 


a(t) = x + f(x) — vjj [x + g{x)\ , (56) 

where f(x) and g(x) are arbitrary C 1 functions of x. The differential 
equation satisfied by x(t) can be obtained from Eq. (30) as 

[x + g(x)\ x - g'{x)x 2 + {f(x) [x + g(x)] - f(x) [g'(x) - l] } x + 

J 2 [x + g{x)} 2 + f 2 (x ) = 0, (57) 

where a prime denotes the derivative with respect to the argument of 
the functions / and g , respectively. Eq. (57) has the first integral 

x = cucot (c ot + a) [x + g(x)\ — f(x). (58) 


3.2.1. Particular case 3: g(x) = (f3 — 1) x, f(x) = —■'yx + 5x n 

As an example of the application of the previous results we assume that 

the functions g(x) and f(x) have the functional forms 

g (x) = (P - 1) x, f(x) = ~jx + 5x n , (59) 


where n,/3,'y,S = constants, and n € N > 1. With these choices 
Eq. (57) takes the form 

f3xx + (1 — /3)x 2 + [{P{n — 1) + 2 )5x n — 2^x\ x — 2y 5x n+1 + 5 2 x 2n + 
x 2 {P 2 uj 2 + y 2 ) = 0, n € N > 1. (60) 

The first integral Eq. (58) becomes 


x = \/3li cot(c vt + a) + 7 ] x — 5x n . 
Therefore the general solution of Eq. (60) is given by 


x{t) = sin /3w {ojt + a) e 7< 
n € N > 1. 


A-5 gin (H-V^iut + a^-Wdt 


(61) 


1/(1 -n) 


(62) 
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3.2.2. Particular case 4- g{x) = 0 
In the case g(x) = 0, Eq. (57) reduces to 

x + — (-^-xf(x) \ x + a?x + ^ — = 0, (63) 

x \dx J x 

Eq. (63) is a standard Lienard type equation. The hrst integral of 
this equation can be obtained from the condition a(t) = a + 
as 

x(t) = u cot(wf + a)x(t) — f(x). (64) 

Equations of the form (63) and (64) were considered in [34], where 
it was shown that Eq. (63) can be exactly solved if the function f(x) 
is of the form [34] 

f{x) = —fix + 5x n , (65) 

with n, (3 and 5 constants. By assuming that 

f{x) = fax 2 + v, (66) 

where g and v are constants, Eq. (63) becomes 

x + + (w 2 + 2/j,v) x + g 2 x 3 + — = 0, (67) 

x x 

with the first integral given by 

x = c jJ cot (cut + a) — /ax 2 — v. ( 68 ) 

Eq. ( 68 ) is a Riccati type differential equation [14]. With the help of 
the transformation 

1 i) , 

x = —, (69) 

H v 

where v is a new C 2 function, Eq. ( 68 ) can be transformed to a linear 
second order differential equation, 

v = UJ cot (c at + a) v — gu. (70) 

By introducing a new independent variable r defined as 

r = — uo cos (cut + a ), (71) 


Eq. (70) takes the form 


d 2 v i-iis 

d? 2 + 1 -t 2 /u> 2V = 1 


(72) 


Therefore the problem of solving the second order non-linear differential 
equation Eq. (67) has been reduced to the problem of solving the linear 
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second order differential equation Eq. (72). The general solution of 
Eq. (72) is given by 


u(r) 


1 1 /I- PVT 1 1 r 


Ci 2 F\ ^ 4/i^w 2 + 1 

^T' <T2jF > (i - Iv'W+i, ^4^2 + l + i; ^ 


+ 

)• 

(73) 


where C\ and C 2 are arbitrary constants of integration. 

It is important to note that for v = 0 Eq. (67) reduces to the 
Modified Emden Equation Eq. (5), whose properties have been studied 
in detail in [24]-[26] and [29, 30]. 


3.2.3. Particular case 5: /(x) = 0 
In the case f(x ) = 0, Eq. (57) becomes 


x — 


9 '( x ) -2 , ,2 


x + ur [x + g(x)\ = 0, 


x + g{x)' 
with the first integral given by 

x = oj cot (ut. + a) [x + g(x)\. 
Therefore the general solution of Eq. (74) is given by 

dx 


exp 


= y4sin(o;f + a). 


x + g{x) J 

In the case g(x) = gox n , go, n constants, Eq. (74) takes the form 

gonx n_1 .2,2, 2 n n 

X - ttX + U X + U) (7nX = 0, 

1 + g 0 x n ~ 2 

with the general solution given by 

x 


(1 + gox' 1 


- —r, —-V = Asin(u;t + a). 

-ixl/fn-i) v > 


(74) 


(75) 


(76) 


(77) 


(78) 


3.3. a(t) = x(t) + f (x) — ic 0 [x(t) + g (x)] 
With the choice 


a(t) = x(t) + f (x) - i uj [x(t) + g (x)], (79) 
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where / (x) and g (x) are arbitrary functions of x, the solution gener¬ 
ating equation Eq. (30) is equivalent to the following non-linear gener¬ 
alized Lienard type equation, 

{ [/' (x) + 1 ]g(x) + x [/' (x) + 1 ] - [/ (x) + x) g' (x)} x + / (x ) 2 + 
xf (x) + w 2 [g (x) + x ] 2 = 0. (80) 

Eq. (80) has the first integral 

x = lo cot (ujt + a) [x + g (x)] — / (x). (81) 


3.3.1. Particular case 6: g(x) = 0 

In the particular case g(x) = 0, we obtain the differential equation 

[fi±)+l]i+ npi + im +w i x=0 . (82) 

XX 

The general solution of Eq. (82) can be obtained by using the first 
integral 

x = ui cot (ut + a)x — / (x). (83) 

As an example we consider that the function / (x) is of the form 

f{x) = ~x + b t (84) 


with b an arbitrary constant. For this choice of /(x), Eq. (82) takes 
the form 


1 b x 3 x 2 b 2 o 

-x- -- - —— + — + w x = 0 , 
4 2 x lb x x 


with the general solution given by 


(85) 




x(t) = — { 26sin [2(u;t + a)] + 86sin 3 (wf + a) coslcot + a) + 
3w I 


3cioj sin 4 (cut + a) >, 


( 86 ) 


where c\ is an arbitrary constant of integration. 


3.3.2. Particular case 7: / (x) = 0 
In the case / (x) = 0, we obtain 

[g (x) + x — xg' (x)] x + J 2 [g (x) + x ] 2 = 0 , 


and 


x = lo cot (oot + a) [x + g (x)], 


(87) 

( 88 ) 
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respectively. With g (.x) = cx, Eq. (87) becomes 

x 2 

x + c 2 co 2 - -h 2coj 2 x + lo 2 x = 0, (89) 

X 

with the general solution given by 

cuj(a-\-tuj) 1 

x(t) = Ae c2 “ 2 +! [| clo cos (ut + a) — sin (ut + a)|] c2 “ 2 +! , (90) 

where A is an arbitrary constant of integration. 


3.4. a(t) = x + f (t, x, x) — iuj [x + g (t, x, .4)] 


Finally, we consider the general deformation of the space phase coordi¬ 
nates x and x by assuming a general dependence of the functions / and 
g on all coordinates t, x and x, respectively. The analysis of this case 
is straightforward, and represents a generalization of all the particular 
cases of deformations previously considered. 

Therefore we summarize the main results of our paper by formulat¬ 
ing the following 

Theorem. Let / = / (t, x, x) and g = g (t, x, x) be two C l arbitrary 
functions of the independent variable t. of the dependent variable x, 
and of its derivative x, respectively. Then the second order ordinary 
differential equation for x, 


df (t,x,x) 
dx 

•v 

dg (t,x,x) 
dx 

dg (t, x, x) 
dt 


+ 1 


>x + 


g (■ t , x, x) + x 
df (t,x,x 


df (t,x,x) 
dx 


dx 

>x + [g (t, x, x) + x\ 


[g(t,,x,x) +x] - f{t,x,x) 


- [f (t, x, x) + x]x 
dg (t, x, x 


dx 


- 1 


(df{t,x,x) 2 . 

< - -Q t -he o [g (■ t , x, X) + x] } - 


f (t, X, X) + f (t , x , i) _ = 0, 


dt, 

has the first integral 


dx 


(91) 


x(t) = uj cot (uit + a) [x(t) + g (t, x, x)] — / (t, x, x). (92) 
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4. Further generalizations 


In the present Section we briefly consider the possibility of extending 
the method presented in this paper to the case of other types of differen¬ 
tial equations. This can be done by generalizing the solution generating 
equation Eq. (30) to other integrable cases. For example, the solution 
generating exactly integrable equation 

X(t) = — 2icoX(t) — ibX 2 (t.) — ib, b = constant ^ co, (93) 


generates, for 


m 


x(t) — i wx(f) 
x(t) + i cox(t) ’ 


the non-linear differential equation 


x(t) + 


b x(t ) 2 
CO x(t) 


co(b — oo)x(t) = 0. 


The general solution of Eq. (93) is given by 


(94) 


(95) 


X(t) = 


co + \\Jb 2 — co 2 tanh yjb 2 — co' 2 (a + t ) 


and therefore Eq. (95) has the first integral 


(96) 


x{t) 


= CO 


y/b 2 — co 2 tanh 

y/b 2 — co 2 (t + a) 

+ i (b — co) 

i \Jb 2 — co 2 tanh | 

(^y/b 2 — co 2 {t + a)^j 

1 + b + co 


x{t). 


(97) 


By deforming the phase space coordinates by introducing the func¬ 
tions / and g, from Eq. (93) several classes of non-linear differential 
equations, as well as their first integrals, can be generated. 

Another possibility of generalizing the previous formalism is to con¬ 
sider that the oscillation frequency of the oscillator is a function of 
time, co = co(t), that is, to start with the generating function 


= x(t) - io >(t)x{t) 
x(t) + i co(t)x(t) 


(98) 


Then the generating equation X(t) + 2ico(t)X(t.) = 0 is equivalent 
with the differential equation 


x{t) 


vjt) 

uj{t) 


x(t.) + co 2 (t)x(t) = 0, 


(99) 
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which has the first integral 

x(t) = oo(t) cot 


oj (t 1 ) dt' + a 


x{t). 


(100) 


The time deformation of the phase space coordinates in Eq. (98), 
x(t) —>• x(t)+g(t), x(t) x(t) + f(t) generates the differential equation 


[. g(t ) + x(t)\ x(t) + 


' -n\ 9 (t)u(t) w(f) ' 


x(t) + 


f(t) - + 2g(t)oj 2 (t) I x(t) + u 2 (t)x 2 (t) - f(t)g(t) - 


u(t) 


f(t)g(t)uj(t) 


u(t) 

which has the first integral 


+ f 2 (t) + 9 2 {t)u 2 (t) + g(t)j(t) — 0, 


x(t) = uj(t) cot 


OJ 


(If) dt' + 


a 


[x{t)+g(t)} - f(t). 


( 101 ) 


( 102 ) 


The systematic study of these types of differential equations will be 
performed in a future study. 

In the case of the linear harmonic oscillator the quantity H = aa + 
represents the conserved energy of the system. For the non-linear gen¬ 
eralization of the harmonic oscillator considered in the present paper 
one can define the energy of the oscillator as 

H = acr = {x + f) + u {x + g) = oj . 9 — 

sm (ojt + a 

The time variation of the energy is given by 

dH _ o 2 (x + g){g- f ) 

7 - c\ . 

dt sin (ojt + a) 

It is interesting to note that the energy is conserved not only for / = 
g = 0, but also for deformations of the phase space coordinates that 
satisfy the condition g = /. 


(103) 


(104) 


5. Applications in Physics and Engineering 

As we have already mentioned, Eq. (67) represents an integrable gen¬ 
eralization of the equation of anharmonic oscillations Eq. (5), which 
arises in a large number of physical problems (see [26] and references 
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therein). From a physical point of view Eq. (67) describes the motion 
of a point particles under the influence of dissipation, described by 
a term of the form (3//x + v/x) x, and an exterior anharmonic force 
lF 2 x+fi 2 x 3 + u 2 /x. By using the obtained results several generalizations 
of the Rayleigh equation [12], 

x + F(x) + x = 0, (105) 

describing non-linear systems with one degree of freedom, admitting 
auto-oscillations, can also be solved exactly. 

In the following we will briefly present some physical and engineering 
applications of the methods introduced in the present paper in two ar¬ 
eas of physics and engineering: travelling waves in reaction-convection- 
diffusion systems, and the approximate solution of the equation describ¬ 
ing large-amplitude non-linear oscillations of a uniform cantilever beam 
carrying an intermediate lumped mass and rotary inertia, respectively. 


5.1. Travelling wave solutions in 

REACTION-CONVECTION-DIFFUSION SYSTEMS 


From a mathematical point of view reaction-convection-diffusion equa¬ 
tions represents a particular class of second order non-linear partial 
differential equations, which play a very important role in many areas of 
physics, chemistry, engineering and biology. A large number of distinct 
phenomena, including heat transfer, combustion, reaction chemistry, 
fluid dynamics, plasma physics, soil-moisture, foam drainage, crystal 
growth, biological population genetics, cellular ecology, tumor growth, 
neurology, synergy etc. can be described mathematically by using re¬ 
action - convection - diffusion equations [39] - [43]. The general class of 
reaction - convection - diffusion equations is given by [39, 42] 


du 

~di 


d_ 

dx 


n , \ du 

D{u) if x 


3v 

+ B(u)- + Q(u ), 


(106) 


where u = u(t, x), and the function D(u) plays the role of the diffusion 
coefficient, B(u) may be viewed as a nonlinear convective flux function, 
and Q(u ) is the reaction function, respectively. Particular classes of 
reaction-convection-diffusion equations are the Fisher, or the logistic 
equation, the Newell-Whitehead, or amplitude equation, the Zeldovich 
equation, and the Nagurno or the bistable equation [42]. 

Reaction-diffusion equations have usually the important property of 
admitting travelling wave solutions [39]-[43]. Travelling wave phenom¬ 
ena have been observed, or produced experimentally, in a large number 
of natural reaction, convection and diffusion processes, and they play 
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an important role, from both theoretical and experimental point of 
view, in many fields of physics, chemistry, biology, and engineering. By 
introducing a phase variable £ = x — Vft, where Vf = constant > 0 is a 
constant wave velocity, and by assuming that u = u(x — Vft), Eq. (106) 
takes the form of a second order non-linear differential equation of the 
form 

d 2 u , , / du\ z .du . . „ 

^ 2 +a(u)^—J +/3(u)—+^(u) = 0, (107) 

where 


a(u) 


d_ 

du 


In D(u),/3(u) 


Vf + B(u) 
D(u) 


>7(«) 


Q(u) 

D(u) ’ 


(108) 


Eq. (107) has the same mathematical form as Eq. (57). Therefore 
it follows that if the reaction function D(u), the convection function 
B(u) and the reaction function Q(u) can be represented in terms of 
two arbitrary function f(u) and g(u) as 


a(u) = 


9\u) 

u + g(u)' 


(109) 


f3(u) 

y(u) 


/(u) (g» - 1) 
u + g{u) 


= u?[u + g{u)\ + 


)_ 

u + g{u)' 


( 110 ) 

( 111 ) 


then the general reaction-convection-diffusion equation Eq. (107) has 
the first integral 


du 


oj cot(cu^ + a) [u + g(u)\ — f(u). 


( 112 ) 


Concentration-dependent diffusion, convection, and reaction coeffi¬ 
cients are used to model mass transfer phenomena in polymer mem¬ 
branes [44], as well as water sorption in a polymer matrix composite 
[45, 46]. Several functional forms for the concentration dependence of 
the diffusion, convection and reaction functions have been proposed, 
and used to model various phenomena, including the very simple form 
of a linear dependence of the diffusion function D(u) = D 0 (1 + kpu ), 
where Dq and ko are constants [44]. 


Lienard_JENG.tex; 30/06/2015; 1:39; p.19 



20 


T. Harko and S. D. Liang 


5.1.1. The case g(u) = ((3 — 1 )u and f(u ) = —ju + 5u 2 

As an application of the above results on the exact integrability of the 

reaction-convection-diffusion equations, we consider the simple case 


g (u) = (/3- 1 )u, 


(113) 


and 

f(u) = —ju + du 2 , (114) 

respectively, where (3, 7 and 5 are constants. Thus we obtain 

a(u) = TA. (115) 

0(«) = -y + 6 (l + |) u, (116) 

7 (u) = u 2 I3u + ^ 7 + 5l ^ 1 (II 7 ) 

giving 

D[u) = D o u {1 ~ 0)/13 , (118) 

B(u) = D o u (1 ~ 0)/0 -^1 + d (l + ^ju - V f , (119) 

Q(u) = DquP-MP oj 2 (3u + + 5u ^ , (120) 

where Dq is an arbitrary constant of integration. 


With this choice of the reaction, convection and diffusion functions, 
Eq. (107) has the first integral 


du 

dt, 


= [f3uj cot(w£ + a) + 7 ] u — du 2 


and has the exact travelling wave solution 

sin^w^ + a)e 7 ^ 


«(0 = 


A — 5 f sim w (u;£ + a)e^d$, 


( 121 ) 


( 122 ) 


Hence we have obtained an exact travelling wave solution of the gen¬ 
eral reaction-convection-diffusion Eq. (107), under the assumptions of 
a specified form of D(u ), B{u) and Q(u). If the coefficients (3 and u 
satisfy the condition (3uj = 1, then the general solution of the reaction- 
convection-diffusion equation (107) with diffusion, convection and re¬ 
action functions (118)-(120) is given by 


«(0 


sin(w£ + a)e 7 ^ 

A — (5e 7 ? [ 7 sin(cj^ + a) — u cos(cu£ + a)] / (q 2 + a; 2 ) 


(123) 
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5.2. Large amplitude free vibrations of a uniform 

CANTILEVER BEAM 

Many engineering structures can be modelled as a slender, flexible 
cantilever beam, which at an intermediate point along its span carries 
a lumped mass with rotary inertia [47]. In the linear theory the vi¬ 
brations of a cantilever beam are usually analyzed by introducing the 
approximation that assumes that the frequency of the free vibrations 
is independent of the amplitude of the motion, that is, that the system 
is isochronous. This approximation is valid only when the amplitude 
of the oscillations is relatively small. However, since beam systems 
are slender and flexible, their flexural vibrations often have relatively 
large amplitudes. Due to their mathematical complexity, such nonlinear 
problems do not allow an exact solution, and therefore approximate 
techniques must be used for their study. In [47], to obtain the approx¬ 
imate solution to the period of oscillation of the cantilever beam the 
single-term harmonic balance method, and the two terms harmonic bal¬ 
ance method was used. Several approximate or semi-analytical math¬ 
ematical methods have been developed for the study of the free large 
amplitude nonlinear oscillations of a slender, inextensible cantilever 
beam, like the homotopy perturbation method [48], modified Lindstedt- 
Poincare methods [49, 50], and the Optimal Homotopy Asymptotic 
Method [51]. In [52] six different analytical methods are applied to 
solve the dynamic model of the large amplitude non-linear oscillation 
equation. 

In the following we adopt the non-linear oscillation model introduced 
in [47]. We consider a beam clamped at the base, and free at the tip. 
The beam caries at an arbitrary intermediate point s = d along its 
span a lumped mass M having a rotary inertia J. To simplify or 
mathematical model we introduce the conservative assumption that 
the beam is uniform, having constant length l, and mass m per unit 
length, respectively. As compared to its length, the thickness of the 
beam is assumed to be small. Moreover, in the following we ignore the 
effects of the rotary inertia and of the shearing deformations. Also the 
beam is assumed to have the important property of being inextensible. 
Note that this implies that the length of the beam neutral axis remains 
constant during the entire oscillatory motion. Then, with the help of 
the above assumptions, the potential energy V. due to the bending of 
the beam, can be written as 


V = 



(124) 
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where El is the modulus of flexural rigidity, £ = s/7 is the dimensionless 
arc length, and R(£,t) is the radius of curvature of the beam neutral 
axis. The kinetic energy T of the beam is given by 



( x 2 + y 2 ) di + -M ( x 2 + y 2 ) + - J 0 


52 


rjl 


(125) 


where y = d/l , a dimensionless quantity, is the relative position param¬ 
eter of the attached inertia element, and 8 is the slope of the elastica, 
given by sinO = dy/ds. In order to simplify the mathematical analysis 
one expands first the kinetic and potential energies into power series, 
with the nonlinear terms retained up to fourth order. Then, by looking 
for an approximate single mode solution of the form ?/(£, t) = $(£)«(£), 
one obtains by simple calculations the discrete single mode, single 
coordinate cantilever beam Lagrangian as 


L = — (oqu 2 + az\ 2 u 2 u 2 — f3f( 


\ 2, ,4\ 


(126) 


where oq, aq, a ,3 , a 4, A, and j3\ are constants. After a rescaling of the 
time coordinate and of the constants i = 1,2,3,4, A and /3, we 
obtain the equation of motion corresponding to the Lagrangian (126) 
as [47], [51] 


u + a- 


u 


• 2 i 2 

;U + U 


U 


+ P- 


1 + au 2 \ 1 + au 2 1 + au 2 

We write the above equation in the form 


u + a 


u 


1 + au 2 


■ 2 , 2 
U + L) 


U + 


u 3 (/? — a) 
1 + au 2 


= 0 . 


= 0, 


(127) 


(128) 


where a = a^/ (p 2 aq) and (5 = 2a 4 / (p 2 a 2 ), where p 2 = Q//3 is the 
dimensionless frequency, and Q is the frequency of the assumed mode 
of the associated linear beam [47]. We try to obtain some approximate 
solutions of Eq. (128) by using the analogy between Eq. (128) and 
Eq. (74). In order to exploit this analogy we first introduce the function 
g(u), satisfying the differential equation 


au d'( u ) 

1 + u 2 u + g' 

with the general solution given by 


(129) 


g(u) 


1 

2 


2yfac\ + sinh 1 (y/au) 
y/a( 1 + au 2 ) 


(130) 
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where c\ is an arbitrary constant of integration. In order to obtain an 
approximate solution of Eq. (128) we tentatively try to approximate 
the function 


h(u) 


((5 — a)u 3 

(1 + era 2 ) ’ 


(131) 


by g(u). To see if such an approximation is valid, we expand both g{u) 
and h(u) in power series near u = 0. We thus obtain 


1 2 era 3 3 n 4 4 a 2 u 5 

q(u ~ ci- ac\u -1 —a c\u H- 

' 2 38 15 


—a 3 citt 6 + O 

16 


h{u ) « (/? — a)it 3 + a (a — /3)u 5 + O (it 7 ) . 
For ci = 0 Eq. (132) becomes 


(« 7 )> 

(132) 

(133) 


9{u) 


au 3 

IT 


+ 


4 a 2 u 5 

15 


+ O (u 7 ) . 


(134) 


The function — g(u ), as given by Eq. (132), is not a good approxima¬ 
tion of h(u) for arbitrary (3. However, for the particular case f3 = 2/3 a 
we obtain 


h{u) F 

s + 0(u 5 ), 

(135) 

g{u) f 

sa ^ +0 (u 5 ). 

(136) 


Therefore for this particular case, and in this order of approximation, 
Eqs. (128) and (74) are the same. Hence in Eq. (136) the ~ sign could 
be substituted by the = sign. 

Thus, by approximating h(u) by g(u), it follows that the solution of 
Eq. (128) can be represented in an integral form as 


A sin(o;f + a). (137) 


Eq. (137) can be used as a starting point for obtaining some power 
series representations of the solutions of the large-amplitude free vibra¬ 
tions of the considered slender inextensible cantilever beam, which is 
assumed to be undergoing planar flexural vibrations. On the other hand 
it is important to mention that Eq. (137) may not hold for arbitrary f3. 

From a physical point of view the terms proportional to a in Eq. (128) 
are of inertia type. They are due to the kinetic energy of the axial 
motion, and they arise as a result of using the approximation of the 
inextensibility condition. The first of these two non-linear terms has a 


exp 


sinh 1 (ySu) 


■ du 


\/a(au 2 + 1 ) 


+ U 
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softening effect, while the second has a hardening effect [47]. On the 
other hand the term proportional to /3 in Eq. (128) is of a hardening 
static type, due to the potential energy stored in bending [47]. There¬ 
fore the condition f3 = (2/3)a can be interpreted as giving a specific 
condition relating the potential energy stored in the bending with the 
kinetic energy of the axial motion. 


6. Discussions and final remarks 

In the present paper we have considered a class of differential equations 
that can be obtained by considering a deformation of the phase space 
coordinates of the linear harmonic oscillator, while keeping unchanged 
the solution generating equation of the harmonic oscillations. The de¬ 
vised method is an application of the quantum mechanical approach to 
the harmonic oscillations to the classical regime. The basic quantities 
used in the present approach are the ’’creation” and ’’annihilation” 
functions a and a + , which in quantum mechanics are represented by 
Hermitic operators [35], and the solution generating equation, con¬ 
structed with the help of these functions. By deforming the phase space 
of the harmonic oscillator by adding two arbitrary functions / and 
g to the position and velocity coordinates, respectively, the solution 
generating equation leads to several classes of dissipative, non-linear 
ordinary second order differential equations, representing a non-linear 
generalization of the harmonic oscillator. The non-linear generaliza¬ 
tions of the linear oscillations can be constructed by using a systematic 
approach based on the adopted form of the ’’deformation” and of the 
’’annihilation” and ’’creation” functions. All the resulting equations 
have a first integral given by the first order differential equation 

X f 

- = u cot (c A + a ). (138) 

x + g 

An interesting question, intrinsically related to the methods devel¬ 
oped in the present paper, is the important question of the isochronous 
nature of the solutions, that is, the independence of the amplitude of 
the oscillation on its frequency. As was shown in [26], all equations 
that can be transformed to the linear harmonic oscillation equation do 
have this property. On the other hand, for all the solutions obtained 
in the present paper the angular frequency c v of the oscillations of the 
’’deformed” system continues to remain the same as that of the linear 
oscillation. Therefore all the obtained solutions exhibits the property 
of amplitude independence of the frequency of oscillation, and the de¬ 
formation of the phase space of the simple linear harmonic oscillator 
generates isochronous oscillating systems. 
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As some simple examples of applications of the obtained results 
in physics and engineering we have briefly discussed the possibility 
of obtaining exact travelling wave solutions of the general reaction- 
convection-diffusion systems, and approximate solutions for the nonlin¬ 
ear equation describing the large amplitude free vibrations of a uniform 
cantilever beam. The travelling wave equation for the general reaction- 
convection-diffusion systems belongs to the type of equations studied in 
the present paper, and as such it can be interpreted mathematically as 
being related to a phase space deformation of the harmonic oscillator 
equation. For variable diffusion, convection and reaction functions, sat¬ 
isfying some given relations, a large number of travelling wave solutions 
can be obtained. We have also presented a particular exact travelling 
wave solution, obtain by prescribing some specific functional forms for 
the functions D, B and Q. 

Most of the mathematical studies of the reaction-convection-diffusion 
equations have been done by using either asymptotic evaluations, or 
qualitative methods. The usual way to study nonlinear reaction-diffusion 
equations similar to Eq. (107) consists of the analysis of the behav¬ 
ior of the system on a phase plane (du/d£,u). This approach, useful 
for qualitative analysis, is insufficient for finding any exact solution, 
or for testing a numerical one. Finding exact analytical solutions of 
the complicated reaction-convection-diffusion models can considerably 
simplify the process of comparison of the theoretical predictions with 
the experimental data, without the need of using numerical procedures. 

Even it presents strong similarities with the differential equations 
studied in the present paper, the equation (128) describing the nonlin¬ 
ear oscillations of uniform cantilever beams does not belong to the class 
of equations that can be obtained by deforming the linear harmonic 
oscillator equation. However, approximate semi-analytical solutions of 
the equation can be obtained under the assumption that the generalized 
anharmonic term, proportional to u 3 , can be approximated by the same 
function g that appears as the coefficient of the term u 2 in the exact 
solution. Of course the validity of such an approach must be checked 
by performing a full comparison between the exact numerical solution, 
and the different order approximations obtained by using the procedure 
introduced in the present paper. 
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